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We study the effects of integrability breaking perturbations on the non-equilibrium evolution 
of many-particle quantum systems. We focus on a class of spinless fermion models with weak 
interactions. We employ equation of motion techniques that can be viewed as generalizations of 
quantum Boltzmann equations. We benchmark our method against time dependent density matrix 
renormalization group computations and find it to be very accurate as long as interactions are weak. 
For small integrability breaking, we observe robust prethermalization plateaux for local observables 
on all accessible time scales. Increasing the strength of the integrability breaking term induces a 
“drift” away from the prethermalization plateaux towards thermal behaviour. We identify a time 
scale characterizing this cross-over. 


In classical mechanics, integrable few-particle systems 
can be understood in terms of periodic, non-ergodic 
motion in action-angle variables. Breaking integrabil¬ 
ity by adding a weak perturbation induces a fascinat¬ 
ing crossover between integrable and chaotic motion, 
which is described by the celebrated KAM theory [1], 
In essence, classical few-particle systems with weak inte¬ 
grability breaking retain aspects of integrable motion on 
intermediate time scales. Recently it has emerged, that 
similar behaviour occurs in the non-equilibrium evolu¬ 
tion of isolated many-particle quantum systems. Start¬ 
ing with the seminal work of Rigol et al [2] it has become 
clear that there is a dramatic difference between the late 
time behaviour of isolated integrable and non-integrable 
quantum many particle systems prepared in initial states 
that are not eigenstates of the Hamiltonian. Generic sys¬ 
tems thermalize [2-13], i.e. exhibit relaxation of local 
observables towards a Gibbs ensemble with an effective 
temperature, while integrable systems evolve towards a 
generalized Gibbs ensemble [2, 12-35]. Starting with the 
work of Moeckel and Kehrein [36] it was then realized 
that models with weak integrability breaking perturba¬ 
tions exhibit transient behaviour, in which local observ¬ 
ables relax towards non-thermal values that retain infor¬ 
mation of the proximate integrable theory. This has been 
termed prethermalization , and has been established to oc¬ 
cur in several models [36-46]. Crucially, it was recently 
observed in experiments on ultra-cold bosonic atoms [47- 
49] . The general expectation is that prethermalization is 
a transient effect, and at “sufficiently late times” non- 
integrable systems thermalize. While this appears nat¬ 
ural, there is scant evidence in support of this scenario. 
The reason is that available numerical [50] or analyti¬ 
cal [36, 41, 45] methods are not able to reach late enough 
times. The exception is the case of infinitely many dimen¬ 
sions, where it was shown in a particular example that 
a weakly non-integrable model thermalizes [51]. Here we 
address these issues in the context of weakly interacting 


one dimensional many-particle systems. This case has 
the important advantage that the accuracy of approxi¬ 
mate methods can be benchmarked by comparisons with 
powerful numerical methods like the time dependent den¬ 
sity matrix renormalization group (t-DMRG) [50]. More¬ 
over, the existence of many strongly interacting one di¬ 
mensional integrable systems makes it possible to verify 
that the qualitative behaviour we find persists for arbi¬ 
trary interaction strengths. 

We focus on the weak interaction regime U < J\ of the 
three-parameter family of spinless fermion Hamiltonians 


L, 

H(J 2 ,S,U) = -J! X [ x + (- 1 )' 5 ] (4 c z+i + H.c.) 
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Here q and cj are spinless fermion operators on site i 
and the hopping amplitudes describe nearest-neighbor 
and next-nearest-neighbor hopping respectively, while 
0 < S < 1 is a dimerization parameter. Finally there 
is a repulsive nearest neighbour interaction of strength 
U. From here onwards we set J\ = 1 and measure all 
the energies in units of J±. There are a several limits in 
which (1) becomes integrable: (i) U = 0 describes a free 
theory; (ii) 5 = J 2 = 0 corresponds to the anisotropic 
spin-1/2 Heisenberg chain [52]; (iii) the low-energy de¬ 
grees of freedom for J 2 = 0 and 5, U <C 1 are described 
by the quantum sine-Gordon model [53]. Away from 
these limits, the model is non-integrable. Our protocol 
for inducing and analyzing non-equilibrium dynamics is 
as follows. We prepare the system in an initial density 
matrix po that is not an eigenstate of H(J 2 , S, U ) for any 
value of U. We then compare the expectation values 
of local operators for time evolution with the integrable 
H(J 2 ,6,0) and (weakly) non-integrable H{J 2 ,5,U) re¬ 
spectively. For (7 = 0 our model is non-interacting, and 
concomitantly in the thermodynamic limit expectation 
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values of local operators relax to time independent val¬ 
ues described by a generalized Gibbs ensemble. In the 
following we analyze how a small integrability breaking 
interaction U > 0 changes the non-equilibrium evolution. 
We stress that our protocol differs in a very important 
way from the weak interaction quenches analyzed previ¬ 
ously [51, 54]. In these works there is no dynamics at 
all for U = 0. Hence quenching the interaction from 
zero to a finite value simultaneously breaks integrability 
and induces a time dependence into the problem. This 
masks the interaction induced modification of the inte¬ 
grate post-quench dynamics. Quantum quenches in the 
model (1) with J 2 = 0 were previously studied in Ref. [41] 
by numerical and analytical methods. On the accessi¬ 
ble time scales robust prethermalization was observed, 
but no evidence for eventual thermalization was found. 
While our manuscript was being completed a paper ap¬ 
peared in which techniques similar to the ones we employ 
here were used to analyze quantum quenches in the case 
8i = 5f = 0 [54]. No prethermalization in our sense 
was observed for the aforementioned reason that there is 
no dynamics without integrability breaking in this case, 
but instead evolution towards a thermal steady state was 
found. Given that U is small, a convenient basis for ana¬ 


lyzing quench dynamics is obtained by diagonalizing the 
quadratic part of the Hamiltonian. This results in 

0) = EE e a (k)al(k)a a (k), (2) 

a=± k >0 

where a± ( k ) are momentum space annihilation op¬ 
erators obeying canonical anticommutation relations 
{a a (k),a'p(q)} = S a ^S k , q , and e a (k) = -2J 2 cos(2fc) + 
2 ay/S 2 + (1 — 5 2 ) cos 2 (fc) are single particle dispersion 
relations of the two bands of fermions. The system is 
initially (at time t = 0) prepared in a density matrix p 0 , 
and subsequently evolves according to 

p(t) = e -iH( J 2,6f,U)t poe iH(J 2 , Sf ,U)t' (3) 

Using equation of motion (EOM) techniques [51, 55] anal¬ 
ogous to the ones employed in derivations of quantum 
Boltzmann equations [56, 57], we obtain evolution equa¬ 
tions for the two-point functions 

n a p(q, t ) = Tr[p(t)al(q)ap(q)}. (4) 

The EOM can be cast in the form 


’^ j otft{]^i ^)— it a. ft T 4zC/e" ^ ^ ^)^7i fti}^'> 0) ^ft^\ 0) 

71 

-U 2 f d t'^2 ^2 K 2p( k i’ k 2-,k-,t-t')n^ 2 (ki,t')n 73 ^(k 2 ,t') 

•'° T fc 1 ,fe 2 > 0 

-U 2 dt'^2 E L ap(k 1 ,k 2 ,k 3 ;h,t-t')n^ 2 (k 1 ,t')nr l3li (k 2 ,t')n~ /5l6 (k 3 ,t'), (5) 

7 fci,fc 2 ,fc3>0 


where e a p(k) = e a (k) — ep(k). Explicit expressions for 
the kernels J, K , L and details of our derivation are 
given in the Supplemental Material. The solution of the 
set of integro-differential equations (5) is numerically de¬ 
manding. We designed an algorithm that scales as L 3 xT 
where T is the number of time steps and L the number of 
lattice sites. This allows us to reach long times Jit ~ 80 
on large systems L ~ 320 (a similar scaling was proposed 
in Ref. [54]). Given the expectation values (4), we may 
readily calculate the single-particle Green’s function 

G{j,l\t) = Tr[p(f)cjcJ 

= ^E E 7a(fcJb/3(MK*/3(M), (6) 

k> 0 a,ft=± 

where the coefficients 7 a {k,j) are given in the Supple¬ 
mentary Material. A crucial check of the accuracy of 
our approach is provided by a direct comparison to pre¬ 
vious t-DMRG computations [41], In Fig. 1 we present 


a comparison of Q(L/2,L/2 + 1) between EOM and t- 
DMRG results for a quench where the system is prepared 
in the ground state of H( 0, 0.8,0) and time evolved sub¬ 
ject to the Hamiltonian 77(0,0.4,0.4). We see that even 
for relatively large U = 0.4, there is excellent agree¬ 
ment between the two methods for all times accessible 
by t-DMRG. Similar levels of agreement are found for 
other Q(L/2, L/2 + j) with j = 2,3,4, 5. This agree¬ 
ment suggests that the EOM method is very accurate for 
small values of U and short and intermediate time scales. 
The advantage of the EOM method is that it allows us 
to access later time scales than the t-DMRG computa¬ 
tions reported in Ref. [41]. As long as the interaction 
strength U is sufficiently small, we observe very long- 
lived prethermalization plateaux, as is exemplified in the 
inset in Fig. 1. There, the thermal value has been com¬ 
puted by quantum Monte Carlo simulations on a system 
with L = 100 sites. 

In order to investigate if and how the prethermalized 
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FIG. 1. (Color online) ^ + 1 ;t) for a quench where the 
system is prepared in the ground state of H (0, 0.8, 0) and time 
evolved with 77(0, 0.4, 0.4) for a system with L = 256 sites. 
The EOM results (red line) are in excellent agreement with 
t-DMRG computations [41] (circles). Inset: prethermalized 
behaviour persists over a large time interval. 



FIG. 2. (Color online) C?(F, — 1; t) for a system with Hamil¬ 

tonian 77(7 2 , 0.1, 0.4) and sizes L = 360, 320 initially prepared 
in a thermal state (7) with density matrix p(2, 0,0,0). The 
expected steady state thermal values are indicated by dotted 
lines, while the black dashed lines are exponential fits to (8). 

regime evolves towards thermal equilibrium it is conve¬ 
nient to invoke a non-zero 7 2 . In essence, 7 2 allows us to 
tune the cross-over time scale between the two regimes. 
In order to access the dynamics for a larger range of en¬ 
ergy densities we consider thermal initial density matrices 
of the form 

g-/3 H{J 2 ,S,U) 

= rj y^ e -/3 H (j 2 ,s tU )y C 7 ) 

Figs. 2 and 3 show results for the time evolution of 
the Green’s function for a system prepared in the ini¬ 
tial state (7) with density matrix p(2, 0,0,0), and time 
evolved with Hamiltonian 77( J 2 , 0.1,0.4). In contrast to 



FIG. 3. (Color online) Real (Inset: imaginary) part of 
C/(F, F + 2; t) for a system with Hamiltonian 77(72,0.1,0.4) 
and sizes L = 360, 320, that was initially prepared in a ther¬ 
mal state (7) with density matrix p(2,0,0,0). The expected 
steady state thermal values are indicated by dotted lines, 
while the black dashed lines are exponential fits to (8). 

the case J 2 = 0, U = 0.4, we now observe a slow drift 
towards a thermal steady state. Increasing J 2 enhances 
the drift. The thermal values shown in Figs. 2 and 3 
are obtained as follows. The energy density is given 
by e = Tr[p(2, 0, 0,0)77(7*2,0.1,0.4)]/L and determines 
the effective temperature 1 / f3 e g of the thermal ensemble 
for the post-quench Hamiltonian 77(72,0.1,0.4) through 
e = Tr[p(/3 eff ,7 2 ,0.1,0.4)77(J 2 ,0.1,0.4)]/L [58], We de¬ 
termine /3 e g by exact diagonalization of small systems up 
to size L = 16, and then use the same method to com¬ 
pute the single-particle Green’s function in thermal equi¬ 
librium at temperature l//3 e g. We note that Q is 
real for odd separations \i — j\. For even \i — j\ the imagi¬ 
nary part is non-zero but small and relaxes towards zero. 
We find that the observed relaxation towards thermal 
values is compatible with exponential decay 

~ Q(i,j) th + A tJ (j 2 ,5, , (8) 

where G(i,j)th is the thermal Green’s function at temper¬ 
ature l/Pefi [59]- The decay times Ty(7 2 ,5, U) are quite 
sensitive to the value of 7 2 . This can be understood by 
noting that large values of 7 2 modify the band structure 
of the non-interacting model by introducing additional 
crossings at a fixed energy. This, in turn, generates ad¬ 
ditional scattering channels that promote relaxation. 

A natural question is whether the integral equation 
(5) can be simplified in the late time regime by removing 
the time integration, in analogy with standard quantum 
Boltzmann equations (QBE) [56, 57]. Here we are faced 
with the difficulty that the structure of our EOM (5) 
is quite different from the ones studied in Refs [56, 57]. 
However, in the case Sf = 0 numerical integration of the 
full EOM (5) suggests that the “off-diagonal” occupation 
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numbers become negligible at late times n.\ _( k , t) ss 0 

and it is possible to derive a QBE for “diagonal” occu¬ 
pation numbers. The QBE for 5f = 0 reads 

n aa (k,T) = -^2 ^2 KZt(p,q'’k)n^{p,T)n ss (q,T) 

7 , S p,q >0 

“E E LZ*(.P,q,r;k) n 77 (p, r)nss{q 1 T)n ee {r, r). (9) 

7,<5,e P,q,r >0 

Here r = C/ 2 i is the usual rescaled time variable, to 
1/U is the time at which the kinetic equation is initial¬ 
ized and the functions K , L are given in the Supple¬ 
mental Material. The QBE agrees with the EOM for 
sufficiently late times (an example is shown in Fig. 4, see 
the discussion below). Because of its simpler structure, 
the QBE allows us to access later times than we able to 
reach with the EOM approach. In particular, employ¬ 
ing the QBE we conclude that for weak interactions the 
relaxation times in (8) scale as [60] 

T- j 1 (J 2 ,5 f =0 : U)cxU 2 . (10) 

This is in contrast to the I/ 4 scaling found for interaction 
quenches in the infinite dimensional Hubbard model [51]. 

To establish more comprehensively that the integra- 
bility breaking perturbation leads to thermalization, we 
consider the (Bogoliubov) mode occupation numbers 
n a / 3 (q,t) themselves. The mode occupation operators 
are not local in space, and hence it is not a priori clear 
that their expectation values should eventually thermal- 
ize; see however Ref. [61]. Importantly, we only con¬ 
sider initial states with finite correlation lengths, which 
implies that Q(j,l;t) are exponentially small in |j — l\ 
as long as \j — l\ Jit [62]. This, together with the 
fact that Q(j,l',t) decay exponentially fast in time for 
\j — l\ < Jit , suggests that n a p(q,t) should relax in 
the regime 1 -C J\t -C L. In Fig. 4 we present the 
mode occupation numbers n aa (k 1 t ) at several different 
times for a system of size L = 320 prepared in the den¬ 
sity matrix p(2, 0,0.5, 0) and evolved with Hamiltonian 
#(0.5,0,0.4). For short and intermediate times J\t < 70 
we use the full EOM, while late times are accessible only 
to the QBE. The QBE is initialized at time to = 20, 
and is seen to be in good agreement with the full EOM 
until the latest times accessible by the latter method. 
We observe that at intermediate times both n+_t_(fc,f) 

and n ( k,t ) slowly approach their respective thermal 

distributions at the effective temperature l//3 e fi' intro¬ 
duced above. The “off-diagonal” occupation numbers 
n + _(k,t), calculated by integrating the full EOM, ap¬ 
proach their thermal value zero in an oscillatory fashion. 
The observed behaviour of the mode occupation num¬ 
bers strongly suggests that the weak integrability break¬ 
ing term indeed induces thermalization. 

We note that in the QBE framework the final relax¬ 
ation is towards the non-interacting Fermi-Dirac distri¬ 
bution with an effective temperature set by the kinetic 




FIG. 4. (Color online) Occupation numbers n++(fc, t) and 

n (k , t ) initialized in the thermal state (7) p( 2,0, 0.5, 0), and 

time evolved with #(0.5, 0, 0.4). The solid lines are the results 
of the EOM (L = 320) for various times. The dotted lines 
are computed by means of the QBE (L = 320). The black 
solid line is the thermal value found by means of second order 
perturbation theory in U. 

energy at the time the Boltzmann is initialized [57, 63], 
signalling the importance of corrections to the QBE at 
very late times. Such corrections, arising from higher 
cumulants, are important for obtaining the power law 
behaviour expected at very late times (for certain observ¬ 
ables) after quenches in non-integrable models [64, 65]. 

In this work we have developed a method that allows 
us to analyze the effects of a weak integrability break¬ 
ing interaction on the time evolution of local observables 
after a quantum quench. We have shown that there is 
a crossover between a prethermalized regime, character¬ 
ized by the proximity of our model to an integrable the¬ 
ory, and a thermal steady state. The observed drift of 
G(i,j',t) in time towards its thermal value is exponential 
and characterized by a time scale proportional to U~ 2 . 
The models considered here feature a global U( 1) symme¬ 
try (particle number conservation). A preliminary anal¬ 
ysis suggests that the scenario found here, a prethermal¬ 
ized regime followed by a cross-over to a thermal steady 
state, occurs also in absence of this C7(l) symmetry [60]. 
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DIAGONALIZING THE NON-INTERACTING HAMILTONIAN 

The non-interacting part of the Hamiltonian (1) 

8) = —J\ ^][1 + S(— 1) ](Cj C l+1 + C^ +1 C Z ) — J2 ^^( C J c l+2 + C l+2 C l )> 


is diagonalized by the canonical transformation 


c l = —ff E E k \ S ) a <*( k ) 

V L k>0a=± 


where the coefficients are given by 

7 ±( 2 j - l,k\6) = e - ik{2j ~ 1) , 7 ±( 2 i,fc| d) = ±e~ ik2j e -*v>UB) > e ~ivUS) = 


— cos k + iS sin k 
\J cos 2 k + 5 2 sin 2 k 


In the new basis we have 


H 0 (J 2 ,S) = EE e a (k)al(k)a a (k ) 

Q!=± fc>0 


(SI) 


(S2) 


(S3) 


where the single-particle dispersions are 

e a {k) = —2 J 2 cos(2fc) + 2aJ\\jH 2 + (1 — 6 2 ) cos 2 (fc). 

Applying the same transformation to the interaction part of the Hamiltonian H lnt = U^2^ =1 cj C;c|, , C;+i gives 

/tint = U EE V ol (k)al 1 (ki)al l2 (k2)a a3 (k 3 )a ai (k4). (S4) 


fc>Q 


Here we have introduced the notations a = (aq, a 2 , a 3 , 0 : 4 ), k = (Aq, fc 2 , k 3 , Aq) and fc > 0 is a shorthand notation for 
ki > 0 for all i = 1,..., 4. The interaction vertex factor can be written in a conveniently antisymmetrized form 

1 


V a {k) = - - S g n ( P ) S gn(Q)Va PiaQl a P2 a Q2 ( k P 1 i k Q 1 > k P2, k Q 2 ) , 


p,Qgs 2 

e i(k 3 -k 4 ) 


K(k) = ^ 2L 6 kt - ka+ka - kt , 


0 


p i(k 3 — fc 4 ) / v 

- 2L (araae^iWe-^W -03046^3 W e -^( tf )J S kl . k2+k3 . ki± ^ 0 , 


(S5) 


where P = (Pi,P 2 ) and Q = ( Qi,Q 2 ) are permutations of (1,2) and (3,4) respectively. 


EQUATIONS OF MOTION 

The equations of motion (5) are derived by following the steps set out in Ref. [56] for deriving quantum Boltzmann 
equations. The starting point are the Heisenberg equations of motion (EOM) for the fermion bilinears fi a p(q,t) = 
a' a (q,t)ap(q,t) . They are of the form 

d 

—h a p(k,t) = i [H, n a p(k, £)] = i [e a (k,S ) - ep(k,6)]h a p(k,t) + iU ’^'^Y"p(k,q)A a (q,t) , (S6) 

a q> 0 

where we have defined A a (q,t) = a\ Xi {qi,t)o) OL2 {q 2 ,t)a^{q^,t)a OLA {q^t), and 

*?) — ^y0,Q!4 &k, q4 Va. 1020:3 o (q) H - $ fi ,a. 3 $k,q 3 Va.\a.20i0t.4 ; ^Q) ^ o, 02 ^ k^q 2 ^ 0 i\^ 0304 (^ 0 ) ^a,Q!i ^k,qi Vfioi20t.30t4 ((?) • 







In the second step we consider the Heisenberg equation of motion for the operator A a (q,t) 


d r - 1 

—A a (q,t) = i H,A a (q,t) = iE a (q)A a (q,t) + iU^^V^p) 

7 p>0 


A-i (p,t),A a (q,t) 


(S7) 


where E a (q) = e ai (qi) + £ 02 ( 92 ) ~ £ 03 ( 93 ) — £ 04 ( 94 )- Integrating (S7) in time and then taking an expectation value 
with respect to our initial density matrix po> we have 


(A a (q,t)) = (A a (q,0))e itE ~W +iU 


dsj2Y, ei(t ~ s)Ea{q)v ^py 

1 p>0 


A-y(p, s),A a (q,s) 


Substituting this back into (S 6 ) leads to an exact integro-differential equation for the mode occupation numbers 
n a p{k,t) = Tr[poh a / 3 (k,t)\, which takes the form 


n a0 {k,t) =i [e a (k, S) - ep(k, 5)] n a/3 (k, t) +iU ET ^Y“ g (k,q) (A a (q, 0)) e ltEa{q) 


a q> 0 


-U 2 dsY, E (A^p,s)A a (q,s))[Y^(k,q)e i ^ E ^V^p)~( a ,q)^( 1 ,p) 

ct,-y q.p >0 


(S8) 


As Wick’s theorem holds for all initial density matrices po we consider, the expectation value (A a (q, 0 )) can be 
expressed in terms of the mode occupation numbers n a0 (k, 0). The eight-point average in (S 8 ) can be decomposed as 

{A^(p,t)A a (q,t)) = f{{n a p(k,t)})+C[(A^(p,t)A a (q,t))] , 


where the first term is the result of applying Wick’s theorem, and (?[•••] denotes terms involving four, six and eight 
particle cumulants (the eight particle cumulant does not contribute because of the antisymmetric structure of (S 8 )). 
In order to turn (S 8 ) into a closed system of integro-differential equations we now assume that the four and six particle 
cumulants can be neglected at all times. This leads to the following system of equations 

h a p{k,t) = ie a0 (k)n a0 (k, t) + AiUe 1 ^ 13(fc) E J yia( k ', t)n 110 (k, 0) - J 011 {k; t)n ail {k, 0) 

71 

-U 2 f di'E E K2p( k i’k2]k;t-t')n^ 2 {ki,t')n 73 ^(k2,t') 

7 fci,fe 2 >0 

-U 2 di'E E kj 2p(k\,k2, fc 3 ; fc;i — i / )n 7 i 72 (fci,i , )n 7 3 74 (fc 2 ,i')n 76 T 6 (fc 3 ,i'), (S9) 

"'° 7 fci,fe 2 ,fc 3 >o 

where 7 = ( 71 ,..., 75 ) and we introduced the functions 


J a0 (k-t)=e ie ^ ( - k)t EE ^0:7273/3 (&? Q .1 Q.1 k)e n ^ 2 ^ 3 ((?, 0 ). 

7273 q >0 

k$ ,/c4 >0 v,v' 




ki , k- 


;;M) = 8 EE^ 


17376 "75 7472 


In- B- k- t) - Hi ^7173^74175^7672 
^-k 1 k 2 k 1 k 2 -k 3 k 1 k 3 k 1 


( a,/3;k;t .), 


12 fc 4 >0 


x k'q( a ^'^',t) = Y2 p {k\q)V a {q)e tE -' { ' k ) t - (7,fe) O (a,q). (S10) 

The occupation numbers at time f = 0 for a system prepared in an initial state with density matrix p(/3, 0, Si, 0) and 
time evolved with Hamiltonian H{J- 2 ,S / ,U) are readily calculated using Wick’s theorem 


n aa [k) = - - - cos (ipk(Sf) - Vk(8i)) tanh(/3e( l 0 ) (fc)/2), 
n a/3(fc) = Isin^fc^/) - <Pk{8i)) tanh(/3e^°' ) (fc)/ 2 ), 

Here the dispersions £a(k) are given by (S3) with J 2 = 0 and (5 = <5*. 


a = ±, 

(Sll) 

3 • 

(S12) 
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QUANTUM BOLTZMANN EQUATION (QBE) FOR 8 f = 0 


The equations of motion obtained in our case are generally quite different from known cases, in which (matrix) 
QBEs can be derived [56, 57, 63]. An exception is the case Sf = 0, where it is possible to obtain a QBE for the 
“diagonal” occupation numbers n aa (q 7 t) as we will now show. Numerical integration of the full EOM, suggests that 
for small interactions [/ « 1 and at sufficiently late times t > f 0 ~ U — 1 the occupation numbers n aa (k,t) depend 

only on the variable r = U 2 t , while n.\ (fc,t) ~ 0. To describe this late time regime it is convenient take the formal 

scaling limit U — > 0, t —» oo keeping r = tU 2 fixed of the EOM (5). In this limit the EOM take the form 


haa,(/c,T) — hn^4i/7 tU ') , n^ lOL (k 1 0) d a ^ 1 (/^,Tf/ )n Q; ^ 1 (/t,0)[]' 

7i 

+ E [ % e iE ^~°)u~ 2 F (k- fc; a) , (S13) 


fc>0‘ 


where we have collected most of the integrand of the u-integral into a single function F(k; k;a) in order to lighten 

notations. The second contribution on the right hand side can be simplified by using our assumptions that n_|_ (k, a) « 

0 and n aa (k,a) are slowly varying functions of a for <r>U. We thus have 


rda 


,iE(k)(r 


(S14) 


The first term vanishes in our scaling limit. We regularize the integral in the second term by replacing E{k) —> 
E(k ) + irj , 77 is small and positive 

nt 

lim / = —--= D(E(k )). (S15) 

u^oJu-! E(k) + ir] y y y J 

The first contribution on the right hand side depends only on the initial mode occupation numbers n a| g(fc,0). In the 
case Sf = 0, the leading contribution at late times is obtained by evaluating the momentum sums by a saddle point 
approximation. This gives 


4* 

U 


T { ^71 a [k] t) n i lQ: (k, 0) J 0:71 ( k;t .) n e* 7 i (M)} 

71 


sin(e + _( 0 )f) 
Ut 3 / 2 


(A a (fc)e ie +-W‘ + c.c.), 


(S16) 


where A a (k) is an amplitude depending on the initial state and the vertex function. The right hand side of (S16) 
vanishes in the scaling limit. Putting everything together, we obtain the following QBE in the scaling limit 

n aa (k,T) = ~y y KZ 1 a 2 (ki,k 2 ; fc)n 7 l 7 l (fci,r)n 7272 (fc 2 ,r) 

7 fci,fc 2 >0 

~y y Aa 1 a 7273 ( fc i 1 fc 2,fc3;fc)n 7 l 7 l (/ci,r)n 7272 (fc 2 ,r)n 7 3 73 (fc3,r). (S17) 

7 ki,k 2 ,k 3 >0 


Here the kernels are given by 

KT(h, k 2 \ q ) = 4 y y («, m 

&3i&4>0 

Li^(k 1 ,k 2 ,k 3 \ q ) = 8y y 

v ki> 0 


(a, m -wy («, m 


X 


2\q ( a ’PW) = Y ap( k \q) v <x{<l) D (E~,(k)) - (7 ,k) O (a,q). 


(S18) 


When implementing the QBE for finite L , the parameter 77 in (S15) must be kept finite (see e.g. [63]). The results 
presented in this paper are for 77 = 0.0005. 





